/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://gitee.com/shenyczz/csharp_csk
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;

namespace CSharpKit.Numerics.Utils
{
    public class KMatrix : ICloneable<KMatrix>, IEquatable<KMatrix>
    {
        public KMatrix()
            : this(0, 0, null) { }

        public KMatrix(int rows, int cols)
            : this(rows, cols, null) { }

        private KMatrix(in KMatrix rhs)
            : this(rhs.Rows, rhs.Cols, rhs.Datas) { }

        private KMatrix(int rows, int cols, IEnumerable<double> datas)
        {
            Rows = rows;
            Cols = cols;
            Datas = null;

            if (datas != null) Datas = datas.ToArray();
            else if (Rows * cols > 0) Datas = new double[cols * rows];
            else Datas = null;
        }

        // 
        public int Rows { get; private set; }
        public int Cols { get; private set; }
        public double[] Datas { get; private set; }

        internal int Size => Cols * Rows;

        public bool IsEmpty => isEmpty();
        public bool IsPhalanx => isPhalanx();
        public bool IsSymmetry => isSymmetry();
        public bool IsSymmetry3 => isSymmetry3();
        public double Determinant => MatrixDeterminant();




        public double this[int row, int col]
        {
            get => Datas[position(row, col)];
            set => Datas[position(row, col)] = value;
        }

        public void Empty()
        {
            Rows = 0;
            Cols = 0;
            Datas = null;
        }

        public void Zero()
        {
            for (var i = 0; i < Size; i++)
            {
                Datas[i] = 0;
            }
        }





        public void Resize(int rows, int cols)
        {
            var mat_old = this.Clone();

            setSize(rows, cols);

            // 恢复数据
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    if (i < mat_old.Rows && j < mat_old.Cols)
                    {
                        this[i, j] = mat_old[i, j];
                    }
                }

            }
        }


        /// <summary>
        /// 生成n阶单位阵
        /// </summary>
        /// <param name="n">单位阵阶数</param>
        /// <returns>true/false</returns>
        public bool GenerateIdentity(int n)
        {
            var result = true;

            try
            {
                setSize(n, n);

                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        if (i == j)
                            Datas[i * n + j] = (1.0);
                        else
                            Datas[i * n + j] = (0.0);
                    }
                }
            }
            catch (System.Exception ex)
            {
#if DEBUG
                Debug.WriteLine(ex);
#endif
                result = false;
            }

            return result;
        }

        /// <summary>
        /// 生成三对称阵
        /// </summary>
        /// <param name="b">n维数组,存储三对称阵中的主对角线元素</param>
        /// <param name="c">n维数组,存储三对称阵中的次对角线元素</param>
        /// <returns>true/false</returns>
        public bool GenerateSymmetry3(double[] b, double[] c)
        {
            var result = true;

            try
            {
                int n = b.Length;

                Empty();
                setSize(n, n);

                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        if (i == j)
                        {
                            Datas[i * n + j] = b[i];
                        }
                        else if (System.Math.Abs(i - j) == 1)
                        {
                            int m = System.Math.Min(i, j);
                            Datas[i * n + j] = c[m];
                        }
                        else
                        {
                            Datas[i * n + j] = (0);
                        }
                    }
                }

            }
            catch (System.Exception)
            {
                result = false;
            }

            return result;
        }

        /// <summary>
        /// 生成托伯利茨(Toeplitz)矩阵 
        /// </summary>
        /// <param name="tu">托伯利茨矩阵上半阵</param>
        /// <param name="tl">托伯利茨矩阵下半阵</param>
        /// <returns>true/false</returns>
        public bool GenerateToeplitz(double[] tu, double[] tl)
        {
            var result = true;

            // if(fabs(tu[0]) < 1.0e-12)	return false;	//tr[0] != 0

            try
            {
                int n = tu.Length;

                Empty();
                setSize(n, n);

                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        if (i <= j)
                            Datas[i * n + j] = tu[j - i];
                        else    //i>j
                            Datas[i * n + j] = tl[i - j];
                    }
                }

            }
            catch (System.Exception)
            {
                result = false;
            }

            return result;
        }

        /// <summary>
        /// 矩阵转置
        /// </summary>
        /// <returns>转置矩阵</returns>
        public KMatrix Transpose()
        {
            var result = default(KMatrix);

            try
            {
                var old = this.Clone();

                var new_rows = Cols;
                var new_cols = Rows;
                setSize(new_rows, new_cols);

                for (var i = 0; i < new_rows; i++)
                {
                    for (var j = 0; j < new_cols; j++)
                    {
                        Datas[position(i, j)] = old[j, i];
                    }
                }

                result = this;
            }
            catch (System.Exception)
            {
                throw;
            }

            return result;
        }












        /// <summary>
        /// 全选主元法求矩阵行列式的值
        /// </summary>
        /// <returns>行列式的值</returns>
        public double MatrixDeterminant()
        {
            // 不是方阵
            if (!isPhalanx())
                return 0;

            var m = this.Clone();

            // 矩阵阶数
            var stRank = Rows;

            // 全选主元
            double MaxValue;
            int iSign = -1, jSign = -1;   // 主元的位置标志
            double Det = 1;      // 行列式的值
            int nSgn = 1;       // 符号
            for (int k = 0; k < stRank - 1; k++)
            {
                MaxValue = 0.0;

                for (int i = k; i < stRank; i++)
                {
                    for (int j = k; j < stRank; j++)
                    {
                        // 求m(i,j)绝对值
                        var tmp = System.Math.Abs(m[i, j]);
                        if (tmp > MaxValue)
                        {
                            iSign = i;			//记下主元位置
                            jSign = j;
                            MaxValue = tmp;
                        }
                    }

                }

                // 绝对值最大元素为0则行列式为0
                if (System.Math.Abs(MaxValue) < 1.0e-12)
                    return 0;

                // 绝对值最大元素不在当前行
                if (iSign != k)
                {
                    nSgn = -nSgn;

                    // 交换行
                    for (int j = k; j < stRank; j++)
                    {
                        swap(m[k, j], m[iSign, j]);
                    }
                }

                // 绝对值最大元素不在当前列
                if (jSign != k)
                {
                    nSgn = -nSgn;

                    // 交换列
                    for (int i = k; i < stRank; i++)
                    {
                        swap(m[i, jSign], m[i, k]);
                    }
                }

                // 对角元相乘
                Det *= m[k, k];

                for (int i = k + 1; i < stRank; i++)
                {
                    if (System.Math.Abs(m[k, k]) < 1.0e-12)
                        throw new Exception("");

                    // 消元因子
                    double d = m[i, k] / m[k, k];

                    // 将主元下方元素消为0
                    for (int j = k + 1; j < stRank; j++)
                    {
                        m[i, j] -= d * m[k, j];
                    }
                }
            }

            // 返回行列式数值 
            return Det * nSgn * m[stRank - 1, stRank - 1];
        }








        #region 操作符重载

        public static KMatrix operator -(KMatrix lhs)
        {
            return lhs * (-1);
        }

        public static KMatrix operator +(KMatrix lhs, double v)
        {
            var mat_new = new KMatrix(lhs);
            var length = lhs.Rows * lhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat_new.Datas[i] = lhs.Datas[i] + v;
            }

            return mat_new;
        }
        public static KMatrix operator +(double v, KMatrix rhs) => rhs + v;

        public static KMatrix operator -(KMatrix lhs, double v)
        {
            var mat_new = new KMatrix(lhs);
            var length = lhs.Rows * lhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat_new.Datas[i] = lhs.Datas[i] - v;
            }

            return mat_new;
        }
        public static KMatrix operator -(double v, KMatrix rhs)
        {
            var mat_new = new KMatrix(rhs);
            var length = rhs.Rows * rhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat_new.Datas[i] = v - rhs.Datas[i];
            }

            return mat_new;
        }

        public static KMatrix operator *(KMatrix lhs, double v)
        {
            var mat_new = new KMatrix(lhs);
            var length = lhs.Rows * lhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat_new.Datas[i] = lhs.Datas[i] * v;
            }

            return mat_new;
        }
        public static KMatrix operator *(double v, KMatrix rhs) => rhs * v;

        public static KMatrix operator /(KMatrix lhs, double v)
        {
            var mat_new = new KMatrix(lhs);
            var length = lhs.Rows * lhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat_new.Datas[i] = lhs.Datas[i] / v;
            }

            return mat_new;
        }
        public static KMatrix operator /(double v, KMatrix rhs)
        {
            var mat_new = new KMatrix(rhs);
            var length = rhs.Rows * rhs.Cols;
            for (int i = 0; i < length; i++)
            {
                if (rhs.Datas[i] == 0) mat_new.Datas[i] = double.NaN;
                else mat_new.Datas[i] = v / rhs.Datas[i];
            }

            return mat_new;
        }

        public static KMatrix operator +(KMatrix lhs, KMatrix rhs)
        {
            //行列相等才能相加
            if (lhs.Rows != rhs.Rows || lhs.Cols != rhs.Cols) throw new Exception("");

            var mat = new KMatrix(lhs);

            var length = lhs.Rows * lhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat.Datas[i] += rhs.Datas[i];
            }

            return mat;
        }
        public static KMatrix operator -(KMatrix lhs, KMatrix rhs)
        {
            //行列相等才能相加
            if (lhs.Rows != rhs.Rows || lhs.Cols != rhs.Cols) throw new Exception("");

            var mat = new KMatrix(lhs);

            var length = lhs.Rows * lhs.Cols;
            for (int i = 0; i < length; i++)
            {
                mat.Datas[i] -= rhs.Datas[i];
            }

            return mat;
        }

        public static bool operator ==(KMatrix lhs, KMatrix rhs)
        {
            var result = lhs.Rows == rhs.Rows && lhs.Cols == rhs.Cols;

            if (result)
            {
                var length = lhs.Rows * lhs.Cols;
                for (int i = 0; i < length; i++)
                {
                    if (lhs.Datas[i] != rhs.Datas[i])
                    {
                        result = false;
                        break;
                    }
                }
            }

            return result;
        }

        public static bool operator !=(KMatrix lhs, KMatrix rhs)
        {
            return !(lhs == rhs);
        }

        public static KMatrix operator *(KMatrix lhs, KMatrix rhs)
        {
            return _MatrixMultiply(lhs, rhs);
        }


        #endregion



        public override bool Equals(object obj)
        {
            return obj is KMatrix matrix &&
                   Rows == matrix.Rows &&
                   Cols == matrix.Cols &&
                   EqualityComparer<double[]>.Default.Equals(Datas, matrix.Datas);
        }

        public override int GetHashCode()
        {
            return HashCode.Combine(Rows, Cols, Datas);
        }

        public bool Equals(KMatrix other)
        {
            return Equals(other as object);
        }

        public KMatrix Clone()
        {
            return new KMatrix(this);
        }



        private int position(int row, int col) => row * Cols + col;

        public void setSize(int rows, int cols)
        {
            Rows = rows;
            Cols = cols;
            Datas = new double[cols * rows];
        }

        private static KMatrix _MatrixMultiply(in KMatrix lhs, in KMatrix rhs)
        {
            // 断定左边矩阵的列数与右边矩阵的行数相等 (左列 == 右行)
            if (lhs.Cols != rhs.Rows) throw new System.Exception("");

            // 生成矩阵新对象，用lhs的行作为新阵的行数，用rhs的列数作为新阵的列数
            var mat = new KMatrix(lhs.Rows, rhs.Cols);

            for (int i = 0; i < lhs.Rows; i++)
            {
                for (int j = 0; j < rhs.Cols; j++)
                {
                    for (int k = 0; k < lhs.Cols; k++)  //左列
                    {
                        mat[i, j] += lhs[i, k] * rhs[k, j];
                    }
                }
            }

            return mat;
        }


        //是否空矩阵
        bool isEmpty()
        {
            return (Rows == 0 || Cols == 0 || Datas == null);
        }

        //是否方阵
        bool isPhalanx()
        {
            return (!isEmpty() && (Rows == Cols));
        }

        //是否对称阵
        bool isSymmetry()
        {
            // 必须是方阵
            if (!isPhalanx()) return false;

            for (int i = 0; i < Rows; i++)
            {
                for (int j = 0; j < Cols; j++)
                {
                    if (System.Math.Abs(this[i, j] - this[j, i]) > 1.0e-12)    // (i,j) == (j,i)
                        return false;
                }
            }

            return true;
        }

        //是否三对称阵
        bool isSymmetry3()
        {
            // 必须是对称阵
            if (!isSymmetry())
                return false;

            // 方阵阶数
            int nRank = Rows;

            for (int i = 0; i < nRank; i++)
            {
                for (int j = 0; j < nRank; j++)
                {
                    if (System.Math.Abs(i - j) > 1)
                    {
                        double dTemp = this[i, j];
                        if (System.Math.Abs(dTemp) > 1.0e-12)  // dTemp!=0
                        {
                            return false;
                        }
                    }
                }
            }

            return true;
        }

        //是否正定阵
        bool isRegular()
        {
            // 正定矩阵必须是方阵
            if (!isPhalanx())
                return false;

            for (int k = 0; k < Rows; k++)
            {
                // 对角元必须大于0
                if (!(this[k, k] > 0))
                    return false;
            }

            for (int k = 2; k < Rows; k++)
            {
                var m = new KMatrix(k, k);
            }

            return true;
        }

        /*
        bool KMatrix::IsRegular() const
        {
            long double ldDet;
            for(k = 2; k <= GetRowNum(); k++)
            {
                KMatrix m(k, k);	//生成一matrix对象

                for(size_t i=0; i<k; i++)
                {
                    for(size_t j=0; j<k; j++)
                    {
                        m(i, j) = GetAt(i, j);	//初始化
                    }
                }

                //ldDet = MatrixDeterminant(m);				// 顺序主子式的值
                ldDet = m.MatrixDeterminant();				// 顺序主子式的值

                if(fabs(ldDet)<1.0e-12 || ldDet < 0.0)	//ldDet <= 0
                    return false;					//不是正定阵
            }
            //---
            return true;
        }
        */


        /*
        ///////////////////////////////////////////////////////////////////////////
        // 名    称: MatrixRank
        // 功    能: 求矩阵的秩: 全选主元高斯(Gauss)消元法
        // 访    问: public
        // 参    数:
        // 返    回: 矩阵的秩
        // 注    释: 
        ///////////////////////////////////////////////////////////////////////////
        size_t KMatrix::MatrixRank()
        {
            size_t i,j,k;

            size_t nRank = 0;						//矩阵秩数

            size_t iSign, jSign;					//主元的位置标志

            size_t stRow = GetRowNum();			//取矩阵行数
            size_t stCol = GetColNum();			//取矩阵列数

            size_t ColRowMin = min(stRow, stCol);	//取行或列最小值

            KMatrix m(*this);				//生成一 KMatrix 对象，用this初始化

            for(k = 0; k < ColRowMin; k ++)	// 全选主元
            {
                long double MaxValue = 0;
                for(i = k; i < stRow; i ++)
                {
                    for(j = k; j < stCol; j ++)
                    {
                        long double tmp(fabs(m(i, j)));	//求m(i,j)绝对值
                        if(tmp > MaxValue)
                        {
                            MaxValue = tmp;
                            iSign = i;					//记下主元位置
                            jSign = j;
                        }
                    }
                }

                //子阵元素绝对值最大者为0
                if(fabs(MaxValue)<1.0e-12) //MaxValue == double(0)
                    break;	//return nRank;
                else				
                    nRank++;		//子阵元素绝对值最大者不为0，矩阵秩加1

                if(k ==(ColRowMin - 1))		//已到最后一行(列)
                    break;	//return nRank;

                if(iSign != k)				//主元不在当前行
                {
                    for(size_t j = k; j < stCol; j ++)	//交换行元素
                        _swap(m(k, j), m(iSign, j));
                }

                if(jSign != k)				//主元不在当前列
                {
                    for(size_t i = k; i < stRow; i ++)	//交换列元素
                        _swap(m(i, jSign), m(i, k));
                }

                for(i = k + 1; i < stRow; i ++)
                {
                    const double d(m(i, k) / m(k, k));		//消元因子
                    for(size_t j = k + 1; j < stCol; j ++)	
                        m(i, j) -= d * m(k, j);		//当前主元右下阵元素作变换
                }
            }
            //-----------------------------------------
            return nRank;
        }


        */








        static void swap<T>(T d1, T d2)
        {
            T dTemp = d1;
            d1 = d2;
            d2 = dTemp;
        }


        // {{END}}
    }

    public class XArray
    {
        // 
        // public 
    }


}
